﻿<?php

//调用案例:
//建国门桥 到 天安门
//echo GPSDistance(39.908551,116.435509,39.908024,116.397743);

//获取角度的弧度制表示
function GetRad($deg){
    return (M_PI*$deg/180);
}


/**
 * 1
 * @param float $lat
 * @param float $lon
 * @param int $type  10最小的经度 11最大的经度 20最小的纬度  21 最大的纬度 
 */
function getMaxLat($lat,$lon,$type=10){
	$offset =0.0005;
	if ($type ==20){
		$lat_to = $lat-$offset;
		$lon_to = $lon;
	}elseif ($type == 21){
		$lat_to = $lat+$offset;
		$lon_to = $lon;
	}elseif ($type == 10){
		$lat_to = $lat;
		$lon_to = $lon-$offset;
	}elseif ($type == 11){
		$lat_to = $lat;
		$lon_to = $lon+$offset;
	}
	

	do{
		$dst = GPSDistance($lat, $lon, $lat_to,$lon_to);
//echo "=".$dst."=<br>";
	
		if ($dst>1000 && $dst<1100)
			break;
		else{
			if ($dst>1100){
				$lon_to-= $offset;
				if ($type ==10)
					$lat_to += $offset;
				elseif ($type==11)
					$lat_to -= $offset;
				elseif ($type==20)
					$lon_to += $offset;
				elseif ($type==21)
					$lon_to -= $offset;
			}elseif ($dst<1000){

				if ($type ==10)
					$lat_to -= $offset;
				elseif ($type==11)
					$lat_to += $offset;
				elseif ($type==20)
					$lon_to -= $offset;
				elseif ($type==21)
					$lon_to += $offset;
			
			}
			
		}
	}while (true);
	
	if ($type ==10 || $type ==11){
		return $lat_to;
	}elseif ($type == 20||$type == 21){
		return $lon_to;
	}
}


//gps 位置计算法
function GPSDistance($lat1, $lon1, $lat2,$lon2) {
    
    $a= 6378137.0;
    $b= 6356752.3142;
    $f= 1 / 298.257223563;
    $rlat1 = GetRad($lat1);
    $rlon1 = GetRad($lon1);
    $rlat2 = GetRad($lat2);
    $rlon2 = GetRad($lon2);
    $L = GetRad($lon2 - $lon1);
    $U1 = atan((1 - $f) * tan($rlat1));
    $U2 = atan((1 - $f) * tan($rlat2));
    $sinU1 = sin($U1);
    $cosU1 = cos($U1);
    $sinU2 = sin($U2);
    $cosU2 = cos($U2);
    $lambda = $L;
    $lambdap = 0.0;
    $iterLimit = 100;
    $sinLambda=0.0;
   	$cosLambda=0.0;
    $sinSigma=0.0;
    $cosSigma=0.0;
    $Sigma=0.0;
    $sinAlpha=0.0;
    $cosSqAlpha=0.0;
    $cos2SigmaM=0.0;
    $C=0;
    do {
        $sinLambda = sin($lambda);
        $cosLambda = cos($lambda);
        $sinSigma = Sqrt(pow(($cosU2 * $sinLambda),2) + pow(($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda), 2));
        If ($sinSigma == 0.0){
            return 0.0;
        }
        $cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
        $Sigma = atan2($sinSigma, $cosSigma);
        $sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
        $cosSqAlpha = 1.0 - pow($sinAlpha,2);
        $cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;
        If (!is_numeric($cos2SigmaM)){
        	$cos2SigmaM = 0;
        }
        $C = $f / 16.0 * $cosSqAlpha * (4.0 + $f * (4.0 - 3.0 * $cosSqAlpha));
        $lambdap = $lambda;
        $lambda = $L + (1 - $C) * $f * $sinAlpha * ($Sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2.0 * pow($cos2SigmaM, 2))));
      	$iterLimit = $iterLimit - 1;
    } While ((abs($lambda - $lambdap) > 0.000000000001) && ($iterLimit > 0));

    If ($iterLimit == 0){
    	return 0.0;
    }
    $uSq = $cosSqAlpha * (pow($a , 2) - pow($b , 2)) / (pow($b , 2));
    $aA = 1 + $uSq / 16384.0 * (4096.0 + $uSq * (-768.0 + $uSq * (320.0 - 175.0 * $uSq)));
    $bB = $uSq / 1024.0 * (256.0 + $uSq * (-128.0 + $uSq * (74.0 - 47.0 * $uSq)));
    $deltaSigma = $bB * $sinSigma * ($cos2SigmaM + $bB / 4.0 * ($cosSigma * (-1 + 2 * pow($cos2SigmaM , 2)) - $bB / 6.0 * $cos2SigmaM * (-3 + 4 * pow($sinSigma , 2)) * (-3 + 4 * pow($cos2SigmaM , 2))));
    $s = $b * $aA * ($Sigma - $deltaSigma);
    
    return $s;
}
//echo GPSDistance(39.908551,116.423509,39.908551,116.435209);
//echo GPSDistance(39.908551,116.435509,39.908551,116.447209);
//echo"<br>";
//echo 116.447209-116.435509;
//echo "<br>";
//echo getMaxLat(39.908551,116.423509);
?>